R Code for GCB-19-0810

rm(list = ls())  # Remove all objects from memory

Influence of most intense ENSO events on the coral cover in the ETP

library(ggplot2)
library(gridExtra)

ENSO_colours<-c("blue", "red") # Type of anomaly
my_colours<-c("#feb24c93", "#fb6a4a93", "#9e9ac893", "#a6d96a93") # Thermal regime

theme_set (theme_classic() + theme(panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(), 
                              axis.line = element_line(colour = "black"),
                              legend.position="none",
                              axis.text.x = element_text(angle = 90, vjust = 0.5),
                              plot.title = element_text(size=12, face="bold"),
                              # panel.border = element_blank(), 
                              panel.border = element_rect(colour = "black", fill=NA, size=1)))

ENSO intensity in region 3

Figure 1a

# Plot ENOS 1971 - 2014
  ENOS<-read.table(file="http://www.cpc.ncep.noaa.gov/data/indices/ersst5.nino.mth.81-10.ascii",header=TRUE)
  ENOS1982_2014 <-subset(ENOS, YR>1968 & YR<2015)   # subset from 1970-2014
  ENOS1982_2014$Date<-paste(ENOS1982_2014$YR, ENOS1982_2014$MON, sep = "-" )
  ENOS1982_2014$Date<-paste(ENOS1982_2014$Date, "15", sep = "-" )
  ENOS1982_2014$Date<-as.Date(ENOS1982_2014$Date, format = "%Y-%m-%d")
  ENOS1982_2014$Anomaly <- factor(ifelse(ENOS1982_2014$ANOM.3 >=0, "Positive", "Negative"))
  
  interp <- approx(ENOS1982_2014$Date, ENOS1982_2014$ANOM.3, n=10000)
  DataENSO <- data.frame(Date=interp$x, Anomaly3=interp$y)
  DataENSO$Anomaly <- factor(ifelse(DataENSO$Anomaly3>=0, "Positive", "Negative"))
      
ENSO<-ggplot(data=DataENSO, aes(x=Date, y=Anomaly3))+
  scale_x_continuous("",expand=c(0, 0))+
  scale_y_continuous("SST anomaly (ºC) in Niño.3",expand=c(0, 0))+
  geom_area(aes(fill=Anomaly), alpha=0.7)+
  scale_fill_manual(values=ENSO_colours) + labs(title="a")+
  geom_hline(yintercept=0, linetype="solid", color="black")+
  geom_line(position = "identity", size=0.3)+
  theme(axis.text.x=element_blank(),
        axis.title.x=element_blank(),
        axis.ticks.x = element_blank())+
  annotate("text", x = 1900, y = 2.0, label = "El Niño \n72-73", size=3)+
  annotate("text", x = 5600, y = 2.0, label = "El Niño \n82-83", size=3)+
  annotate("text", x = 11100, y = 2.0, label = "El Niño \n97-98", size=3)
ENSO

# ENSOb<-ggplot(data=ENOS1982_2014, aes(x=Date, y=ANOM.3))+
#   geom_hline(yintercept=0, linetype="solid", color="black")+
#   scale_x_date("", limit=c(as.Date("1969-01-15"), as.Date("2014-12-31")),
#                breaks= "12 months",
#                      expand=c(0, 0))+
#   scale_y_continuous("SST anomaly (ºC) in Niño.3",expand=c(0, 0))+
#   geom_area(aes(fill=Anomaly))+
#   geom_line(position = "identity")+
#   scale_fill_manual(values=ENSO_colours) + labs(title="(a)")
# ENSOb

DHW and coral cover change

# Load packages 
library(tidyverse)
library(dplyr)
library(data.table)
library(RColorBrewer)

Yearly maximum DHW experienced by the best-represented coral reefs sites from 1982 to 2014

  • Select and transform DHW data
# Select DHW data file: maxDHW.csv
  dhw <- read.csv("Data/maxDHW.csv", header=TRUE)

# Exclude years 1981, and 2017-2019
  dhw <- dhw[!dhw$Year == 1981, ]

# Create groups for plot
  ts.Cano_island <- subset(dhw, dhw$Site == "Cano_island")     
  ts.Uva_island <- subset(dhw, dhw$Site == "Uva_Island")          
  ts.Gorgona_island <- subset(dhw, dhw$Site == "Gorgona_Island")          
  ts.Darwin_island <- subset(dhw, dhw$Site == "Darwin_North_Anchorage")          
  ts.Bahia_Banderas <- subset(dhw, dhw$Site == "Bahia_Banderas")

FIGURE 2a

  • Plot max DHW in each location
# Year, maxDHW
  par(oma = c(1,1,1,1)) 
  par(mar = c(4,4,1,1), cex.axis = 1,  las = 1)  
  par(lwd = 1)   
  
  plot(ts.Cano_island$Year,ts.Cano_island$maxDHW,   # Create a white false line with real data
       xlab="Years", ylab="Degree heating weeks (°C-weeks)",  cex=1.5, cex.lab=0.6, col="white",axes = FALSE,
       ylim = c(0, 28), lwd = 1.5) 
  par(new=T)
  rect(1982, 0, 1983.5, 28, col="#fee8c8", border ="#fee8c8")           # Strong ENSO years
  par(new=T)
  rect(1997, 0, 1998.5, 28, col="#fee8c8", border ="#fee8c8")
  par(new=T)
  plot(ts.Cano_island$Year, ts.Cano_island$maxDHW,type="l",              # Thermally stable
  xlab="", ylab="", lwd = 1.5, col="#8c2d04",axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Uva_island$Year, ts.Uva_island$maxDHW, type="l",                # Thermally stable
       xlab="", ylab="", lwd = 1.5, col="#feb24c",axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Gorgona_island$Year,ts.Gorgona_island$maxDHW, type="l",         # Tropical upwelling
       xlab="", ylab="",lwd = 1.5, col="#fb6a4a", axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Darwin_island$Year,ts.Darwin_island$maxDHW, type="l",           # Galapagos
       xlab="", ylab="", lwd = 1.5, col="#9e9ac8",axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Bahia_Banderas$Year, ts.Bahia_Banderas$maxDHW, type="l",        # Seasonal
       xlab="", ylab="",lwd = 1.5,col="#a6d96a", axes = FALSE,
       ylim = c(0, 28))
  y <- c(0,4,8,12,16,20,24,28)                                            # Axis
  axis(2,lwd = 0.5, at = y,cex.axis=0.6)
  box(lwd = 0.5)
  x <-c(1982:2014)
  YearsR = as.character(c(1982:2014))
  axis(1,  at = x, labels = YearsR, lwd = 0.4, cex.axis=0.5, las = 2) 
  abline(h=4, col="#525252", lty=2)
  abline(h=8, col="#525252", lty=2)

# Legend
  sites <- c("Thermally stable - Caño Island","Thermally stable - Uva Island",
             "Tropical upwelling - Gorgona Island","Equatorial Upwelling - Darwin Island", 
             "Seasonal - Bahía Banderas")
  par(xpd=TRUE)     
  legend(x = 2000, y = 20, 
         xpd = NA,  # or TRUE 
         inset=c(-0.5,0),
         cex = 0.4,
         legend = sites,         # box categories 
         col = c("#8c2d04","#feb24c","#fb6a4a","#9e9ac8","#a6d96a"),
         lty= 1,
         lwd = 1.5,
         x.intersp = 0.5,
         y.intersp = 2, 
         bty="n")

 dev.off() # reset parameters for next plot 
## null device 
##           1

Relationship of coral cover rate of change and thermal stress

  • The DHW index measures the accumulated thermal stress above a bleaching threshold for 12 consecutive weeks.
  • At each site, the maximum monthly mean sea surface # temperature (MMM) is calculated from a long-term satellite climatology dataset.
  • Any temperature surpassing the MMM by 1ºC or more is added for 12 consecutive weeks,
  • because temperatures ~1ºC above the MMM are known to cause coral bleaching (Liu et al., 2006).
  • For example, if the water temperature during four consecutive weeks exceeds the MMM in 1.3°C, # 1.1°C, 1.2°C, and 1.4°C, then,
  • the site accumulates 5°C-weeks when these anomalies are summed.

We tested the hypothesis that the strength of the Log_annual_rate_of_change (response variable) is a function of the maxDHW experienced. We had measured the rate of change from four Thermal regimes, and fitted Thermal regimes as a random intercept, and estimate a common slope (change in the rate of change) for maxDHW across all Thermal regimes by fitting it as a fixed effect.

## Load data
  roc.cover <- read.csv("Data/RoC_DHW.csv", header=TRUE)  ### Select file:  RoC_DWH.csv

# Load packages
  library(lme4)
  library(lmerTest)
  library(MuMIn)
  • Model 6: coral cover change and DHW
model6 <- lmerTest::lmer(Log_annual_rate_of_change ~ maxDHW + (1|Thermal_regime), data=roc.cover)
step(model6)
## Backward reduced random-effect table:
## 
##                      Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)
## <none>                             4 -17.738 43.476                     
## (1 | Thermal_regime)          0    3 -27.772 61.545 20.069  1  7.469e-06
##                         
## <none>                  
## (1 | Thermal_regime) ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##        Eliminated  Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## maxDHW          0 0.65054 0.65054     1 128.59   9.873 0.002082 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## Log_annual_rate_of_change ~ maxDHW + (1 | Thermal_regime)
summary(model6) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Log_annual_rate_of_change ~ maxDHW + (1 | Thermal_regime)
##    Data: roc.cover
## 
## REML criterion at convergence: 35.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7962 -0.1472  0.1167  0.4152  1.8621 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Thermal_regime (Intercept) 0.02501  0.1581  
##  Residual                   0.06589  0.2567  
## Number of obs: 133, groups:  Thermal_regime, 4
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)  -0.085763   0.084256   3.155381  -1.018  0.38032   
## maxDHW       -0.016316   0.005193 128.594887  -3.142  0.00208 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr)
## maxDHW -0.165
r.squaredGLMM(model6)   #  returns the marginal and the conditional R2
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help
## page.
##             R2m       R2c
## [1,] 0.05249512 0.3131505
# marginal R2  describes the proportion of variance explained by the fixed factor(s) alone:
# conditional R2 describes the proportion of variance explained by both the fixed and random factors
# Nakagawa, S. & Schielzeth, H. (2013).

FIGURE 2B: Relationship of coral cover rate of change and thermal stress

# Load packages
  library(tidyverse)
  library(dplyr)
  library(data.table)
  library(RColorBrewer)
 
# Create groups for plot
  Seasonal          = subset(roc.cover, roc.cover$Thermal_regime == "Seasonal")          
  Thermally_Stable  = subset(roc.cover, roc.cover$Thermal_regime == "Thermally_Stable")  
  Tropical_Upwelling = subset(roc.cover, roc.cover$Thermal_regime == "Upwelling") 
  Galapagos         = subset(roc.cover, roc.cover$Thermal_regime == "Galapagos_Islands") 

# Plot
  par(oma = c(1,1,1,1)) 
  par(mar = c(4,4,1,1), cex.axis = 1,  las = 1) 
  par(lwd = 1)
  plot (Seasonal$maxDHW, Seasonal$Log_annual_rate_of_change,
        pch = 21, lwd = 0.5, bg = "#a6d96a93", col="black", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1),
        xlab = "Degree heating weeks (°C-weeks)",
        ylab = "Annual rate of change in coral cover")
  par(new=T); 
  plot (Thermally_Stable$maxDHW, Thermally_Stable$Log_annual_rate_of_change,
        xlab="", ylab="", pch = 21, lwd = 0.5, col="black", bg = "#feb24c93", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1))
  par(new=T); 
  plot (Tropical_Upwelling$maxDHW, Tropical_Upwelling$Log_annual_rate_of_change,
        xlab="", ylab="", pch = 21, lwd = 0.5, bg = "#fb6a4a93",col="black", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1))
  par(new=T); 
  plot (Galapagos$maxDHW, Galapagos$Log_annual_rate_of_change,
        xlab="", ylab="", pch = 21, lwd = 0.5, bg = "#9e9ac893",col="black", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1))
  box(lwd = 0.5)
  x <-c(0,5,10,15,20,25,30)
  axis(1, lwd = 0.5, at = x, cex.axis = 0.5)
  y <- c(-1.5,-1.0,-0.5,0,0.5,1)
  axis(2,lwd = 0.5, at = y, cex.axis = 0.5)
  abline(v=4, col="#525252", lty=2); abline(v=8, col="#525252", lty=2); abline(h=0, col="#525252", lty=2); 
  # Axis text I, II, III, IV
  text(3.5, 0.3, "II", cex = 0.6); text(5, 0.3, "I", cex = 0.6)
  text(3.5, - 0.3, "III", cex = 0.6); text(5, -0.3, "IV", cex = 0.6)
  
  par(xpd=TRUE) 
  legend(x = 10, y = -0.5, 
         xpd = NA,   
         inset=c(-0.5,0),
         legend = c("Equatorial Upwelling", "Thermally stable", "Tropical upwelling","Seasonal"),
         pch=21,
         pt.bg = c("#9e9ac8","#fb6a4a","#fb6a4a","#a6d96a"),
         x.intersp = 0.5,
         y.intersp = 2, 
         bty="n")

 dev.off() # reset parameters for next plot 
## null device 
##           1

Descriptives statistics heat stress

No stress

no.stress <- subset(roc.cover, maxDHW < 4)                      #subset of DHW < 4
length(no.stress$maxDHW)/length(roc.cover$maxDHW)               # ratio no stress: total              
## [1] 0.8496241
1- length(no.stress$maxDHW)/length(roc.cover$maxDHW)               # ratio stress: total              
## [1] 0.1503759
positive.rate <- subset(no.stress, Log_annual_rate_of_change > 0)   #subset of Rate > 0
length(positive.rate$Log_annual_rate_of_change)                     # Number of observations in 
## [1] 48
negative.rate <- subset(no.stress, Log_annual_rate_of_change < 0)   #subset of Rate > 0
length(negative.rate$Log_annual_rate_of_change)                     # Number of observations in 
## [1] 58

No stress seasonal

no.stress.seasonal <- subset(no.stress, Thermal_regime == "Seasonal")   
positive.rate <- subset(no.stress.seasonal, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.seasonal, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)   
## [1] 0.07142857

No stress upwelling

no.stress.upwelling <- subset(no.stress, Thermal_regime == "Upwelling")
positive.rate <- subset(no.stress.upwelling, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.upwelling, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)  
## [1] 1

No stress thermally stable

no.stress.Thermally_Stable <- subset(no.stress, Thermal_regime == "Thermally_Stable")   
positive.rate <- subset(no.stress.Thermally_Stable, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.Thermally_Stable, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)  
## [1] 1.2

No stress Galapagos

no.stress.Galapagos_Islands <- subset(no.stress, Thermal_regime == "Galapagos_Islands") 
positive.rate <- subset(no.stress.Galapagos_Islands, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.Galapagos_Islands, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)  
## [1] 1

Supplementary figures

library(pals)
library(dplyr)
library(reshape2)
library(RColorBrewer)

FIGURE S2: a) Number of regions surveyed per year, b): Number of countries surveyed per year

# Barplot number of surveys per year
  obs_num = as.data.frame(table(cover$Year,cover$Region))    
  colnames(obs_num) = c("Year", "Region", "Freq")
  wide_obs_num = dcast(obs_num, Region~Year, value = "Freq")  
## Using Freq as value column: use value.var to override.
  df = data.matrix(wide_obs_num[,2:45])   
  as.table(df)     
##   1970 1971 1972 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
## A    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## B    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## C    0    0    1    0    0    0    0    0    1    1    0    0    1    0
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    1    0    1    0    4
## F    0    0    0    2    3    2    0    0    0    0    0    1    0    1
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    3   10    2   13    4    4    3    5    2    5    2    3    5    3
## I    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## J    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## K    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## L    0    0    0    0    0    0    0    2    0    0    0    0    0    0
## M    0    0    0    0    0    0    0    0    0    0    0    1    0    0
## N    2    6    0    2    1    0    0    0    0    0    0    1    4    2
## O    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## P    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998
## A    0    0    0    0    0    0    0    0    0    1    0    0    0    0
## B    0    0    5    0    0    0    0    1    0    0    0    0    0    0
## C    0    0    1    2    4    0    0    1    0    0    2    2    0    4
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E   12    1    2    0    1    0    0   11    0    2    0    6    0    1
## F    0    1    1    0    1    0    0    0    2    0    0    0    0    0
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    4    5    2    3    1    1    0    4    1    4    4    0    1    1
## I    0    0    0    0    0    0    2    0    0    0    0    0   16   10
## J    0    0    1    1    1    0    3    1    1    0    0    0    2    2
## K    0    0    0    0    0    0    0    0    0    0    0    0    2    1
## L    0    0    0    0    0    0    0    0    0    1    1    4    4    2
## M    0    0    0    0    0    0    0    1    0    0    0    0    0    0
## N    0    1    5    1    2    0    1    1    1    1    1    1    1    1
## O    0    0    0    0    0    0    0    0    1    1    0    0    0    0
## P    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
## A    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## B    0    0    0    1    0    0    0    1    0    0    4    0    0    0
## C    7    4    3    4    6    4    3    6    3    2    3    2    2    3
## D    0    0    0    0    0    0    0    0    0    0    0    4    0    1
## E    1    0    2    0    4    2    2    0    1    0    2    0    0    0
## F    0    0    1    1    0    1    3    6    2    2    0    0    1    0
## G    0    0    0    0    0    0    0    0    0    0    1    0    0    0
## H    0    1    4   14   12    1    1    1    1    0    1    2    0    0
## I    0    0    0    6    3    0    1    1    4    0   17    6    2    0
## J    0    0    0    1    2    0    0    0    0    5    7    0    0    1
## K    1    5    5    3    2    0    0    0    0    0    2    0    0    0
## L    0    0    2    0    1    0    1    3    2    0   17    3    1    1
## M    0    0    0    0    0    1    2   14    8    0    0    0    0    1
## N    1    1    1    1    1    1    0    0    0    0    1    0    0    0
## O    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## P    0    0    0    0    0    2    0    0    0    0    0    0    0    0
##   2013 2014
## A    0    0
## B    1    1
## C    2    1
## D    4    2
## E    0    0
## F    0    0
## G    0    0
## H    0    0
## I    2    2
## J    0    0
## K    0    0
## L    1    3
## M    0    0
## N    0    0
## O    0    0
## P    0    0
  df[ df>1 ] <- 1 
  levels(cover$Region) # How many regions? = how many colours
##  [1] "Clipperton"                 "Cocos_Islands"             
##  [3] "Colombia"                   "Costa_Rica_Central"        
##  [5] "Costa_Rica_Southern"        "Eastern_Galapagos_Islands" 
##  [7] "El_Salvador"                "Gulf_of_Chiriqui"          
##  [9] "Mexican_Central"            "Mexican_Northern"          
## [11] "Mexican_Southern"           "Nicaragua_Papagayo_Zone"   
## [13] "Northern_Galapagos_Islands" "Panama_Bight"              
## [15] "Revillagigedo_Islands"      "Western_Galapagos_Islands"
  Regions <- c(levels(obs_num$Region))
  color.regions = warmcool(17)
  
  # Figure 2a
  par(mfcol=c(2,1), cex.main = 1.5, mar = c(4,4,1,1), cex.lab = 1.5, cex.axis = 1,  las = 1)
  barplot(df,
          ylim=c(0,15),
          col = color.regions,   
          ylab="Number of regions surveyed per year",  
          xlab="Year",
          #lwd = 1,          
          las = 2,           
          axis.lty = 1, 
          cex.names = 0.4,   
          cex.axis = 0.8)
# Legend as a box
  legend( #"bottomright",    
     x = 2, y = 16, 
     xpd = NA,  # or TRUE 
     legend = Regions,          
     fill = color.regions[1:17],   
     inset=c(-0.5,0),
     cex = 0.35,
     text.font=0.5,            
     x.intersp = -0.5,
     y.intersp = 1,
     lwd = 0.1,
     bty="n")

# Figure 2b
obs_num = as.data.frame(table(cover$Year,cover$Country))    # Table per country
colnames(obs_num) = c("Year", "Country", "Freq")
wide_obs_num = dcast(obs_num, Country~Year, value = "Freq")   
## Using Freq as value column: use value.var to override.
dfc = data.matrix(wide_obs_num[,2:45])   
as.table(dfc)       
##   1970 1971 1972 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
## A    0    0    1    0    0    0    0    0    1    1    0    0    1    0
## B    0    0    2    9    0    0    0    2    0    2    0    1    0    4
## C    0    0    0    2    3    2    0    0    0    0    0    2    0    1
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## F    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    5   16    0    6    5    4    3    5    2    4    2    4    9    5
##   1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998
## A    0    0    1    2    4    0    0    1    0    0    2    2    0    4
## B   12    1    7    0    1    0    0   15    0    6    4   10    4    3
## C    0    1    1    0    1    0    0    1    2    0    0    0    0    0
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    1    0    0    0    0
## F    0    0    1    1    1    0    5    1    2    1    0    0   20   13
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    4    6    7    4    3    1    1    2    2    2    2    1    2    2
##   1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
## A    7    4    3    4    6    4    3    6    3    2    3    2    2    3
## B    1    0    7    1    5    2    3    4    3    0   12    7    1    2
## C    0    0    1    1    0    4    5   20   10    2    0    0    1    1
## D    0    0    0    0    0    0    0    0    0    0    1    0    0    0
## E    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## F    1    5    5   10    7    0    1    1    4    5   26    6    2    1
## G    0    0    0    0    0    0    0    0    0    0   11    0    0    0
## H    1    2    2   15   13    2    1    1    1    0    2    2    0    0
##   2013 2014
## A    2    1
## B    6    6
## C    0    0
## D    0    0
## E    0    0
## F    2    2
## G    0    0
## H    0    0
dfc[ dfc>1 ] <- 1 # replace all values greater than one with one
color.countries <- brewer.pal(9, "Oranges")
countries = c(levels(obs_num$Country))
barplot(dfc,
        col = color.countries,   
        ylim=c(0,10),
        ylab="Number of countries surveyed per year",  
        xlab="Year",
        las = 2,           
        axis.lty = 1,      
        cex.names = 0.4,   
        cex.axis = 0.8)
# Legend as a box
legend( #"bottomright",    
   x = 2, y = 10, 
   xpd = NA,  # or TRUE 
   legend = countries,         
   fill = color.countries[1:9],   
   inset=c(-0.5,0),
   cex = 0.35,
   text.font=1,            
   x.intersp = -0.5,
   y.intersp = 1,
   lwd = 0.1,
   bty="n")

 dev.off() # reset parameters for next plot 
## null device 
##           1

FIGURE S3: Variation in the live coral cover for 1970–2014 per thermal regime

Temporal variation in the live coral cover for 1970–2014. * (a) Tropical upwelling regime. * (b) Seasonal regime. * (c) Thermally stable regime. * (d) Equatorial upwelling regime. The smooth line represents the best fit to the time series with a 95% confidence level interval and a span = 0.3.

library(ggplot2)
library(lemon)

aggr.region$Thermal_regime<-factor(aggr.region$Thermal_regime, 
                                  levels=c("Tropical upwelling", "Seasonal", 
                                           "Thermally stable", "Equatorial upwelling"),
                                  labels = c ("a", "b", "c", "d"))
                                   
summary(aggr.region$Thermal_regime)
##  a  b  c  d 
## 75 28 74 25
FigureS3<-ggplot(aggr.region, aes(x=Year, y=Coral_cover)) + 
  theme_bw() + theme(panel.border = element_blank(),   # set ggplot to white background
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(), 
                              axis.line = element_line(),
                              strip.background = element_blank(),
                              strip.text = element_text(face="bold",hjust = 0, vjust = -1))+
  geom_point()+ geom_smooth (span = 0.3) +
  scale_x_continuous(limits=c(1970, 2014), 
                     breaks = seq(1970, 2014, by=5),
                     expand = c(0.01, 0.01), "Year") +
  scale_y_continuous(limits=c(-19, 91), expand = c(0, 0), "Live coral cover (%)") +
  facet_wrap(Thermal_regime~., ncol=1, scales = "free_x", strip.position ="top")
FigureS3

#ggsave(file="FigureS3.svg", plot=FigureS3, width=6.69, height=6.69)